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Smoothed Particle Hydrodynamics 3.1. Introduction 

3.1 Introduction 

The study of astrophysical phenomena presents a multitude of obstacles to the po- 
tential student. In addition to the usual obstacles of understanding the physical 
properties of the system in question, the sheer scale of astrophysical events renders 
laboratory experiments impossible in the vast majority of cases. To this end, it has 
been necessary to assemble a new, numerical laboratory in the form of computational 
simulations, and conduct experiments and analyses via this medium. The growth 
of computing power over the past 80 years, from the Colossus of Bletchley Park's 
Enigma code-cracking efforts in the 1940s, through ENIAC and Los Alamos Na- 
tional Laboratory's modelling of thermonuclear detonations in the 1950s, up to the 
supercomputers of today, has in turn allowed the computational domain to become 
a mainstay of astrophysical experimentation. 

Two principal approaches to computational simulations have evolved to enable 
these numerical simulations. Eulerian methods use geometric grids, either fixed 
or adaptive (the so-called AMR or Adaptive Mesh Refinement codes), with the 
fluid parameters evaluated over the grid cells. Such codes formed the basis of the 
revolution in Computational Fluid Dynamics (CFD) that started in the late 1960s 
and early 70s, and as such they remain the most widely used approach. Applications 
of such codes cover a huge range, from industrial aerodynamics in the automotive 
and aerospace sectors, to stress calculations and solid mechanics for civil engineering 
and architecture, to chemical reaction modelling and protein folding in biomolecular 
models. 

Lagrangian methods on the other hand dispense with fixed points in space and 
instead evolve the fluid equations in a co-moving frame. A common approach is 
to use discrete particles that are carried with the flow - hydrodynamic (and other) 
properties are then evaluated at the particle positions, and are calculated from a 
weighted average of the values on other local particles. In this manner each particle 
is essentially "smoothed" over a finite volume of fixed mass, and in this way these so- 
called Smoothed Particle Hydrodynamics or SPH codes are naturally adaptive with 
density. Although SPH was originally developed by the astrophysics community, it 
too has found uses and applications in a much wider range of fields. In engineering it 
has been applied to dam breaks and atomised oil lubrication flows, while a number 
of physics engines in computer games use SPH as a basis. The community has grown 
to the point where there is now a Europe-wide network of users called SPHERIC - 
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the SPH European Research Interest Communitjlj. This aims to share advances in 
code development across the user community, and to prevent the re-invention of the 
wheel when it comes the solution of known problems. 

Each of these approaches has advantages and disadvantages with respect to the 
other. Generally speaking, AMR codes have a higher resolution for a given number 
of grid cells than an SPH code with an equal number of particles. Furthermore, they 
can be made to adapt to any flow parameter (although this is not always trivial!), 
while SPH adapts primarily with density only. On the other hand, SPH naturally 
handles vacuum boundary conditions, whilst large grids are required with AMR 
codes to prevent the flow disappearing from the edge of the computational domain. 
As SPH is a Lagrangian method, advection of flow properties is inherent, whereas 
this presents problems for AMR codes, and which usually entails an unphysical 
increase in entropy. In a similar manner, SPH codes can be implemented in such 
a manner that they are inherently conservative of mass, momentum and energy, 
and similarly, unless it is explicitly added in shocks, they likewise conserve entropy. 
Nonetheless it is emphatically not true to say that either SPH or grid code methods 
are "better" than the other, simply that the more appropriate approach should be 
chosen for any given problem, and indeed greater confldence in the results will ensue 
if the two methods concur. 

Having said that, throughout this chapter 1 shall however consider only the SPH 
approach, as it is this one that 1 have used to generate all the results discussed 
henceforth. Furthermore, as all the problems 1 have considered have been fully 
three-dimensional, throughout this chapter 1 shall consider only the derivation and 
discussion of SPH in 3D. This chapter is therefore structured as follows: In Sec- 



tion 13.21 1 shall introduce the basic concepts of the SPH method, then in Section 13.31 
this is used to re- write the fluid equations in a manner that can be solved numerically. 
Section 13.41 discusses the dissipative processes required for the correct implementa- 
tion of artiflcial viscosity and the introduction of entropy in shock waves, and then 
in Section 13.51 1 discuss how the SPH formalism may be made more adaptive still by 
the self-consistent inclusion of variable smoothing lengths. As many astrophysical 
problems are strongly influenced by gravitational forces 1 detail how these may be 
implemented in Section 13.61 In Section 13.71 1 briefly summarise the methods used 
to flnd the nearest neighbours, and then in Section 13.81 1 consider how the code is 



evolved forward in time, and various time-stepping criteria. Finally in Section I3T9] 1 
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briefly outline the properties of the code I have used, point the reader in the direc- 
tion of some standard numerical tests used for code evaluation and consider further 
extensions to the method. 

3.2 SPH Basics 

In the following section I shall discuss the derivation of the SPH formalism from 
first principles, showing how a continuous field can be mapped on to (and thus 
approximated by) a series of discrete particles, and the errors involved in this ap- 
proximation. I then show how derivatives may be calculated, and discuss ways in 
which the particles may be suitably smoothed to represent the field. 

3.2.1 Discrete Approximations to a Continuous Field 

We start from the (mathematically) trivial identity 

/(r) = y"/(r')5(r-r')dr', (3.1) 

V 

where /(r) is any (scalar) function defined on a three-dimensional co-ordinate system 
r ranging over a volume V . Similarly, 5{v) is the Dirac delta function, and r' is a 
dummy variable also ranging over V . 

We may generalise the delta function to a so-called smoothing kernel W with a 
characteristic width h (known as the smoothing length) such that 

\imW{Y,h) = 5{Y), (3.2) 

subject to the normalisation 

fw{r,h)dr' = l. (3.3) 

V 

By expanding W{r — r', h) as a Taylor series, it can be shown that for symmetric 
kernels W{r — r', h) = W{r' — r, h), equation 13.11 becomes 

/(r) = J /(r') W{r - r', h) dr' + 0{h'), (3.4) 

V 
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the second order accuracy arising from the vanishing; of t he kernel gradient at r' = r 
(see for instance |Pricdl2005l : lBenall990l : lMonaghanlll992l ) . Note that more elaborate 



kernels accurate to 0{h'^) can be constructed, b ut these suffer from the problem that 



m 



W{r, h) can become negative in certain ranges (iPricd . l2005l : iMonaghanl . Il992l ) . thus 
potentially leading to negative density evaluations in certain pathological situations. 
Nonetheless for a second order, symmetric kernel, for any finite density p(r) 
within V, equation 13.41 is exactly equivalent to 



/« 



V 



m 

p(r') 



W{r-r',h) p{r')dr' + 0{h^ 



(3.5) 



Discretising this continuous field on to a series of particles of (potentially variable) 



mass m = p(r')dr', the original identity equation 13. II becomes 



/(r) 



J2—f{r.)W{r-r,,h), 

. Pi 



(3.6) 



where now /(r,), rrii and pi = p{vi) are the scalar value, mass and density of the 
jth ^qtj-iiqIq^ an(j I ranges over all particles within the smoothing kernel. Equation 
13.61 therefore represents the discrete approximation to the continuous scalar field 
/ at position r in the computational domain V , and is thus the basis of all SPH 
formalisms. Note that the position r at which the function / is approximated is 
completely general and is not restricted to particle positions, although in practice 
this is where the values are actually evaluated. 



3.2.2 Spatial Derivatives and Vector Calculus 

In order for the SPH discretisation of a field to be useful as a method of solving fluid 
flows, it is clear that the spatial derivatives of any given quantity must also have a 
suitable approximate fornu- Here therefore, I summarise the SPH approximations 
for various vector calculus quantities. 



3.2.2.1 Gradient of a Scalar Field 

The approximation for the gradient of a scalar field can be derived by taking the 



spatial derivative of equation 13.11 and applying the smoothing kernel. Noting that 



^Temporal derivatives are naturally also required, and these will be discussed in due course 
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V = d/dr, we therefore see that 

V/(r) = |-|/(r')5(r-r')dr' (3.7) 

V 

= -^ [ f(r')W(r-r',h)di' + 0(h'^), (3.8) 

or J 



in a similar manner to equation 13.41 Given that the only part to depend on r is the 
smoothing kernel W, and again introducing the density p(r') in both the numerator 
and the denominator we obtain 

V/(r) = I j^-^Wir - r', h)p{r') dr' + Oih'). (3.9) 

V 

Finally this may be discretised in the same way as before, to give 

V/(r) ^ Y. — /(^^) ^^(^ - ^- ^) (3.10) 

as an estimator for the gradient of a scalar field /(r). Notable from the above result 
is that the gradient of a scalar field can be approximated by the values of the field 
itself along with the gradient of the kernel. Computationally this is very useful as 
at no point does V/ have to be evaluated for any particle, whilst the gradient of 
the kernel will be known explicitly for any sensible choice of W. 

3.2.2.2 Divergence of a Vector Field 



Although equation 13.11 was given only for a scalar field, a similar identity may be 
given for a vector field F(R), namely 



F(r) = 

V 



/"F(r')(5(r-r')dr', (3.11) 



Taking the divergence of this with respect to r, and noting once again that the only 
term to depend on r is the smoothing kernel we find that the integral approximation 
becomes 

V ■ F(r) = f F(r') ■ \/W{r - r', h)dr' + ^(/i^), (3.12) 
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and thus as before this can be discretised to obtain the approximation 

V ■ F(r) ^ y — F(r,) ■ VWiv - r^, h). (3.13) 

3.2.2.3 Curl of a Vector Field 

By a precisely similar argument, it is possible to show that the curl of a vector F, 
V X F can be approximated using 

V X F(r) ^ y — F(ri) x WWir - r^, /i), (3.14) 

although this is relatively little used unless magnetohydrodynamic (MHD) effects 
are being taken into account. 



3.2.3 Errors 

The approximations given in equations 13. 6[ I3.10[ 13. 131 and 13. 141 encompass both the 
0{h'^) errors of considering only the integral term, and also the errors inherent in the 
discretisation (which arise due to incomplete sampling of the smoothing kernel). In 
the former case we see that the 0{h'^) errors are reduced by decreasing the smoothing 
length, while the discretisation (sampling) errors are minimised by increasing the 
number of particle s within the smoothing kernel. Barring numerical stability issues 



(JRead et al.l . 120101 ) . this discrete approximation is therefore at its most accurate with 
large numbers of particles contained within a small smoothing length. However, this 
must be balanced against the need for computational speed and efficiency, and hence 
there is a compromise to be struck. 

These errors are neatly illustrated by considering the approximations to a con- 
stant function /(r) = 1 and the zero function, which can be obtained by noting that 
with this definition of /, V/(r) = 0. The SPH approximations for one and zero 
therefore become 

1 ^ Y,— W{r-r,,h), (3.15) 



^ ^!i!i vvr(r-ri,/i). (3.16) 

i P' 

Since in neither case does the equation reduce to an identity, we see that there are 
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inherent errors in estimating even constant functions. Nonetheless, with suitable 
choices for the number of particles within the smoothing kernel and the smoothing 
length, these may be kept to an acceptable level. F or a rn o re det a iled derivat i on an d 
discu s sion o f the se errors, th e read er is directed to IPricd (120051 ) ; iMonaghanI ( 1l992[ ) ; 
Benzl Jl99o[ ) and beadet aP J201ok 



3.2.4 Improved Approximations for Spatial Gradients 

Although the approximations given in equations 13. 6[ I3.10[ 13.131 and 13.141 are those 
that arise most readily from the SPH approximation, it is possible to construct other 
estimators for the gradient of a scalar field. For instance, by noting that for any 
quantity /(r) = l./(r), we see that 

V/(r.l) = l.V/(r) + /(r)Vl (3.17) 

and therefore that 

V/(r) = V/(r) - /(r)Vl. (3.18) 

Clearly, since VI = these forms should be identical. From equation 13. 161 however. 
we see that the SPH approximation for VI is non-zero, and thus using equation 13. 181 
we may define another estimate for V/(r) as 

V/(r) = V— /(r.)Vl^(r-r„/.)-/(r)V^VW^(r-r„/i) (3.19) 
i P^ i P^ 

= ^^.ZMlLMvW^(r-r„/i). (3.20) 

i P' 

This approximation clearly has the advantage that it vanishes identically for constant 
functions. 

A more general class of interpolants arises from considering the vector calculus 
identity 

V(/p")=n/p"-iVp + p"V/, (3.21) 

valid for all n G M. This in turn leads to the following identity for V/ 

W = ^[V(/p")-n/p"-Vp]. (3.22) 

Substituting p and /p" into equation 13. 10[ we obtain a general interpolant for V/, 
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such that 

V/(r) = -^Y,m,, (/(r,)p(r,)"-i - nf{v)p{vY~') VW{v - r„ h). (3.23) 

Two instances of this general case turn out to be particularly useful, namely where 
n = 1 and n = — 1. For the former case we obtain 

V/(r) = JL ^ m, (/(r) - /(r,)) VW{r - r„ h). (3.24) 

This is very similar in form to that given in equation I3.20[ with the exception that 
knowledge of the density at r is required a priori. Although no longer anti-symmetric 
in /(r) and /j, it is nonetheless exact for constant functions. 
In the case where n = — 1 we obtain 

V/(r) = p(r)J]m,(-M+ /^^ WW{r-r,,h). (3.25) 

While this form is no longer exact for constant functions, it is commonly used as an 
estimator for the pressure gradient (VP)/p, as it is pairwise symmetric and as such 
ensures conservation of momentum. This is also the form of the gradient that arises 
naturally from a Lagrangian formulation of the fluid equations, as I shall show in 
Section 13.3.21 



3.2.5 Improved Divergence Estimates 

In a similar manner to the gradient, improved estimates can be made for the diver- 
gence of a vector field. By noting that F(r) = l.F(r), the estimate 

V ■ F(r) ^ V — (F(r,) - F(r)) ■ VW(r - n, h). (3.26) 

can be arrived at, which again becomes exact for constant functions. In a similar 
manner to the expansion given in equation 13. 2H a general class of estimates can be 
arrived at by considering the identity 

V ■ (p"F) = p"V ■ F + rip"-iF ■ Vp, (3.27) 
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the n = 1, — 1 cases of which are given by 

V ■ F(r) ^ ^ V m, (F(r,) - F(r)) ■ VW{r - r„ h) (3.28) 

p(r) ^ 



and 



V ■ F(r) ^ p(r) Y: m. {^, + ^^f ■ VW{r -r^,h) 



(3.29) 



respectively. Once again these estimates have the advantages of being exact for 
constant functions in the former case and pairwise symmetric in (V ■ F)/p in the 
latter case. 

3.2.6 Smoothing Kernels 

From the above it is clear that the choice of smoothing kernel is an important one. 
It must by definition obey the criteria set out in equations 13.21 and 13.31 in that it 
must tend to a ^-function as /i — )■ and it must be normalised so the area under the 
curve is unity. For the purposes of calculating the gradients of quantities it is also 
clear that it should have a continuous and well defined first derivative, and from a 
symmetry argument it should be spherically symmetric, and thus depend only on 
r = |r — r'l and h. 

One of the first choices for the smoothing kernel was the Gaussian function, such 
that 

Wir,h) = j^^e-\ (3.30) 

where x = r/h. However, this has the drawback that W^ > for all r, and thus 
all particles within the computational domain contribute. The computational cost 
of such a kernel therefore scales as (9(A^^), where N is the number of particles in 
the simulation. Given that (for purely hydrodynamical quantities) long range forces 
are negligible, it makes sense to restrict the kernel to those with compact support, 
i.e. make them subject to the condition that W{r, h) = where r/h > k for some 
constant k. This means that the computational cost scales as 0{NN^cigh), where 
-^ncigh is the average number of particles within a sphere of radius r = kh about any 
one particle. 



For this reason, cubic spline kernels are often used (see iMonaghan fc Lattanzio 
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19851 for instance), where the kernel is defined as 



W{r, h) 



1 x^ + -x^ < X < 1; 



Trh"^ 







1 < X < 2; 
a; > 2, 



(3.31) 



where x = r/h as in equation 13.301 Here only particles within 2h of the central 
particle contribute to the smoothing kernel, which is spherically symmetric and 

smoothly different iable for all r. A lthough many other kernels are possible (see 

for example) this 



Monaghan 



1992 



Fulk fc Quinn 



1996:lPrice 



2005; 



Read et al. 



2010 



is a commonly used kernel, and is the one present in the code I have used throughout. 
Note from the above the gradient of the kernel is well defined for all values of x, 
such that 



VW{r, h) 



^W{r,h) 

9 n 



-X 



3x < X < 1 



Tlh^ 



4 



(2- 



X] 







1 <x< 2 
X >2 



(3.32) 
(3.33) 



Finally it is worth noting that in general the form of the kernel makes little 
overall difference to the computational speed of the code. This is because most 
codes tabulate the values of both the kernel and its gradient rather than compute 
them directly, and thus the form of the kernel may be as simple or as complex as 
required, even (theoretically at least) to the extent of being non-analytic functions. 

3.3 Fluid Equations 

Given that the SPH formalism has now been put on a sound mathematical footing, 
in this section I shall use it to obtain approximations to the equations governing 
fluid motion, such that they can be used to construct a viable numerical algorithm 
for solving fluid flows. For the dual purposes of brevity and simplicity I shall here 
consider only the case of an inviscid compressible flow in the absence of body forces, 
although the inclusion of both gravity and (artificial) viscosity will be discussed 
in due course. First however, it is useful to summarise the principal equations of 
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motion in their standard conservative form. 

The continuity (conservation of mass) equation is given by 

f + V ■ (pv) = (3.34) 

where as normal, p is the density, t is time and v is velocity. 

The Euler equation gives the equations of motion in the case of an inviscid fluid, 

and encapsulates the conservation of momentum. In the absence of external (body) 

forces it becomes 

dp\ 



dt 



V-(pv®v) + VP = 0, (3.35) 



where P is the fluid pressure and (8) represents the outer or tensor productp. For 
compressible flows it is also necessary to take into account the energy equation, and 
as such the conservation of energy is embodied in the following equation; 

- + V ■ [(« + P)v] = 0, (3.36) 

where u is the specific internal energy and w = |v| is the magnitude of the velocity 
vector. Finally it is worth noting that these five equations (there are 3 components 
to the momentum equation) contain six unknowns (p, the three components of 
velocity v^, Vy and Vz, P and u). Therefore in order to solve the system we require 
a further constraint; an equation of state is required. All the analysis I shall present 
henceforth uses the ideal gas equation of state, where 

P = k(s)p^, (3.37) 

= (7-l)wp, (3.38) 

where 7 is the adiabatic index (the ratio of specific heats), which throughout has 
been set to 5/3, and k,{s) is the adiabat, itself a function of the specific entropy s. 
In the case of isentropic flows, s and (thus k) remains constant. 

I shall now discuss the SPH formulation of each of the continuity, momentum 
and energy equations in turn. Note that again for the purposes of brevity I assume 
that the smoothing length is held constant (i.e. h = 0, where the dot denotes the 
derivative with respect to time), and is equal for all particles. Individual, variable 
smoothing lengths will be discussed in due course. Furthermore, I assume through- 



^The outer product of two vectors may be summarised as A ® B = AB-^ ~ AiBj (in indicial 



notation). 
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out that the mass of each particle is held constant, such that rrii = const, and again 
that all particles are of equal mass. Although it is possible to have individually 
varying particle masses, the code I use does not have this feature, and therefore I 
have not included a discussion of it here. Finally note that from here onwards, all 
the approximations are evaluated at specific particle positions, as this is how the 
SPH algorithm is implemented within particle-based codes. 

3.3.1 Conservation of Mass 



Using equation 13. 6| we see that in the case of the density, the SPH approximation 
becomes very simple, namely that at particle j the density pj becomes 

Pj = ^miW{rj -r^,h), 

' (3.39) 

i 

where we write Wji = W{rj — rj,/i), and where by symmetry, Wji = Wij. Note 
that here and henceforth, as the SPH formalism is a discrete approximation to the 
underlying continuous medium, we assume equality between the estimator on the 
RHS and the SPH quantity on the LHS. 



(3.40) 



Taking; the full time derivative of equation 13.391 we obtain 


dpj ^^ 

dt ^ 


'dWji dvj dWj, dYi dWjidh' 
dvj dt dvi dt dh dt 


and noting that 





dr,- drj dh 

^ = "- d^="- d^ = °' 
we find that the time derivative of density becomes 

^ = J2m, (v, ■ V,W,, + V, ■ V.W,, 

i 



(3.41) 



where we use \ji = Vj — Vj, and where we note that the gradient of the kernel is 
antisymmetric, i.e. that 

ViW.i = -V,Wj,. (3.42) 
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From equation I3.28| we note that the RHS of equation 13.411 is simply an estimator 
of —pjVj • Vj. Hence equation 13.411 becomes 

^ = -p,V,-v,, (3.43) 

which is simply a reformulation of the continuity equation equation 13.341 using the 
Lagrangian time derivative 

^4 + (v.V), (3.44) 

in which the second term accounts for the advection of flow properties through 
the fluid. Therefore we see that the SPH estimate for density equation 13.391 is 
automatically conservative of mass (as long as equation 13.281 is used as an estimate 
for the divergence of velocity) . 

3.3.2 Conservation of Momentum 

Although there are various ways of deriving the equations of motion consistently with 
the SPH framework, a particularly appealing one is to use the Lagrangian formalism. 
As long as the discrete Lagrangian functional preserves the fundamental symmetries 
of the underlying continuous one, this confers the inherent advantages that the 
resulting SPH equations of motion will automatically fulfil the requisite conservation 
laws (through Noether's Theorem) and also that the only approximations made are 
in the discretisation of the Lagrangian itself. 

3.3.2.1 Linear Momentum 

Defined as the total kinetic energy of the system minus the total internal energy (for 
purely hydrodynamical flows), the Lagrangian functional C for the fluid is 



£(r, v) = / -pv ■ \ — pu dr, (3.45) 



where as before, u is the specific internal energy. For later simplicity, we note 
that through the equation of state (equation I3.38P the specific internal energy is a 
function of density and pressure u = u{p, P), which in turn are functions of position. 
This gives u = u{r). Now if we again make the discretisation mi = pdr, the SPH 



14 



Smoothed Particle Hydrodynamics 



3.3. Fluid Equations 



estimate of the Lagrangian becomes 



UATi 



(3.46) 



where i ranges over all particles. 

The equations of motion for particle j are obtained from the Lagrangian through 
the Euler-Lagrange equations, as follows; 



d /dC 



dt \ d\j J dvj 



dC 



0. 



(3.47) 



By considering each of the terms in this equation it is therefore possible to obtain 
an SPH approximation to the equations of motion that remains fully conservative. 
If we therefore consider the derivative of the Lagrangian with respect to the velocity 
at particle j, we find 



dC d ^ fl \ 

t; — = t; — > ?Tij -Vj ■ Vj — Ui(rj) 



m^Vj-, 



(3.48) 



noting that since the velocities are independent the differential is zero unless i = j. 
Considering now the second term in the Euler-Lagrange equation 13.471 we find 
that 



9£ 

dVn 



E 



dui dPi dui dpi 



dPi dvj dpi dvj 



(3.49) 



where we have used equation 13. 38l to obtain the full derivative of the internal energy. 
In the isentropic (dissipationless) case we see that k,{s) is constant, and thus the 
pressure is a function of density only, leading to 






dui dPi dui 



dpi 



dPi dpi dpi\ dvj 

From the equation of state equation 13.381 we find that 

dui dPi dui Pi 
dPi dpi dpi pf ' 



(3.50) 



(3.51) 



and thus the derivative of the Lagrangian with respect to the position of particle j 
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becomes 





y^ Pi dpi 


|3.39|wetindthat 




dpi 




= 





(3.52) 



(3.53) 
(3.54) 



where we take rj^ = |rjfc|, and use the fact that the kernel is spherically symmetric. 
By direct differentiation, 

OT -u 

-^ = i^ij - ^kj)rik, (3.55) 

with fjfc = Yik/vik the unit vector in the direction of Fj^. Substituting this back into 
equation 13.541 we find that 

-rrr = Z^^k^r—iSij - 6kj)r,k (3.56) 

j u '-J' ik 

= 5Z"^fcV,W/,fc (<5., -4j), (3.57) 

k 

where in the second case we have used the fact that d/dvj = Vj. 
With reference to equation 13.521 we find therefore that 

dC P 

— - = - ^ m~ ^ TTikS/jWik {Sij - 6kj) (3.58) 

"^J i Pi k 



P. ^ ^ P 



-mj 



o k i Pi 

i V Pj Pt / 

where we have changed the summation index in the first term to i and used the fact 
that the gradient of the kernel is antisymmetric, i.e. that VjWkj = —^jWjk- Finally, 
by substituting equations 13.481 and 13.601 into equation 13.471 and dividing through by 
the common factor nij, we find that the SPH equations of motion become 



dt 
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Since this equation is pairwise symmetric in i,j, it is clear that the pressure force 
on particle j due to particle i is equal and opposite (due to the antisymmetry of 
the kernel gradient) to the force on particle i from particle j. In this manner, it is 
clear that this formulation of the equation of motion conserves linear momentum by 
construction. 

3.3.2.2 Angular Momentum 

To check that angular momentum L = r x m\ is conserved, we note that its deriva- 
tive with respect to time should be zero. By using equation 13.611 we see that the 
time derivative of the angular momentum of particle j is given by 

dL,- dv,- . , 

— -^ = m,v,- X V,- + m,ri x — -^ (3.62) 



-mj 



S'"^(^ + ^)'"^'''^^^^^' ^^-^^^ 



since by definition Vj x v^ = 0. The total time derivative of the angular momentum 
is therefore given by the sum over all particles j, such that 

^ = - E E ^^-^^ (7 + f ) ^^- ^ ^^^^r (3-64) 



J ^ 



Hence we see that by reversing the summation indices the entire sum is 
antzsymmetric in i and j, i.e. 



— = -^^m,mJ-f + -|jr,x V,W^,,-, (3.65) 

i j V/^i Pj / 

which can only be the case where the total sum is zero. Hence the angular momen- 
tum is constant with time, and thus angular momentum is explicitly conserved. 

3.3.3 Conservation of Energy 

In the case of a purely hydrodynamical flow, the total energy E = pu + pv^/2 is 
given by the sum of the kinetic and internal energies, such that the SPH estimator 
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becomes 

E = ^mi(-Vi-\i + Uij . (3.67) 

i ^ ' 

Clearly, where energy is conserved the time derivative of the total energy should be 
zero. Taking the time derivative therefore, we find that 

dE sr^ ( dvi dui&Pi duidpi\ 



dt ^ ' \ ' dt dPi dt dp, dt 

where we have again used the fact that in the dissipationless case P = P{p) and we 
can therefore amalgamate the latter two terms of the RHS of equation 13.681 using 
equation l3.51l Using also the equation of motion derived above in equation 13.611 and 
dp/dt from the continuity equation 13.411 we therefore find that 

dE v^ v-^ fPi P\^^^^ 

+ E"^5 E^^(^*-^^)-^^-^i^ (3-70) 

i f^i j 

= ^^m,m,Y^v, + ^v,yV,iy,., (3.71) 

where we have again used the fact that the kernel is antisymmetric to obtain equa- 
tion 13.711 Now using the same argument we used to show that angular momentum 
is conserved, we note that equation 13.711 is antisymmetric under a reversal of i and 
j, and thus must be equal to zero. Hence we find that 

^^0. (3.72) 

and therefore that the total energy is also explicitly conserved. 

A corollary of this is that the time derivative of the internal energy is given by 
the second term on the RHS of equation I3.69| such that 

^ = ^^, (3.73) 

= 4 ^ mi (vj - Vi) ■ ViWji, (3.74) 
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and indeed this is how the internal energy is evolved within SPH codes. 

It is worth noting that the formulation of SPH outlined above is therefore explic- 
itly conservative of mass, momentum (in both forms) and energy. Hence, while there 
are inevitably errors inherent in the SPH discretisation of a continuous medium, 
these are the only errors that appear, at least in the case of a dissipationless hydro- 
dynamical flow. 

3.4 Dissipative Effects 

So far we have assumed the fluid flow to be barotropic (i.e. P = P{p)), and poly- 
tropic, with the polytropic index set equal to the adiabatic index 7, the ratio of spe- 
cific heats. This in turn means that the flow is isentropic, and therefore completely 
dissipationless. While this is an adequate approximation for many incompressible, 
inviscid and unshocked compressible flows, it presents serious problems when it 
comes to modelling transonic and supersonic flow regimes, as the conversion of me- 
chanical (kinetic) energy into heat (internal) energy is not correctly captured. The 
problem occurs because at a shock front, flow properties such as the velocity, pres- 
sure, density and entropy change very rapidly, on the order of the mean free path 
of the gas particles. On large scales therefore these changes appear discontinuous, 
and flow solvers that do not resolve the mean free path (which is all of them) break 
down due to the apparently singular flow gradients. 

There are two principal workarounds that allow numerical codes to solve shocked 
flows. One is to use a Riemann solver in a Gudonov-type code (see for instance 



Inutsuka 



2002 



Cha fc Whitworth 



20031 ). but I shall not go into any detail here as 
this is not the approach used in the code I have used. The alternative approach, 
used in the majority of SPH codes, is to broaden the shock across a small number of 
smoothing lengths. This ensures that the flow gradients do not become infinite, and 
gives the correct asymptotic behaviour away from the shock. This latter method is 
implemented by including an artificial dissipative term in the momentum and energy 
equations that is triggered only in the presence of shocks, and it is this method that 
I shall consider here. 
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3.4.1 Standard Artificial Viscosity Prescription 

Due to the fact that by construction, shock capturing through a viscous process is 
an artificial one, there is considerable latitude in the way in which such an artificial 
viscos ity may be implemented. This being said, it must obey the following general 



rules f lvon Neumann fc Richtmyer 



1950 



Rosswog|,l2009|): 



- The flow equations should contain no discontinuities; 

- The shock front should be of the order of a few times the smoothing length; 

- The artificial viscosity should reduce to zero away from the shock front; 

- The Rankine-Hugoniot equations should hold over length scales larger than 
that over which the shock is smoothed, i.e. 



Po + 



Povo 
Povl 



PlVl, 

Pl + 



Pivl 



— + uo + — 
Po 2 



Pi I'l 

— + ^1 + 17, 

Pi 2 



(3.75) 
(3.76) 

(3.77) 



where the subscripts and 1 refer to pre- and post-shock regions respectively. 

The overall conservation of momentum and energy should not be adversely 
affected, while the entropy should rise from the pre- to post-shock regions. 



By considering the SPH approximation to the momentum equation 13.611 where the 
force is based on pairwise addition of terms of the form P/p^, on dimensional grounds 
it seems sensible to consider an artificial viscosity term 11 of the form 



Hoc 



P 



(3.78) 



for some suitable velocity scale v. Ivon Neumann fc Richtmverl fll950l ) suggested a 
viscous term dependent on the squared velocity divergence (which gives an indication 
of the local expansion or contraction of the fluid), which translates into SPH form 
as 

(n.,)NR = ^'^^^'\^-^-\\ (3.79) 



Ptj 



where h represents a characteristic length scale (in SPH this is equivalent to the 
smoothing length), pij is the average density of particles i and j and /3sph is a 
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constant term of order unity. Noting that to first order 



IV ■ V,: 



;/l2 



(3.80) 
(3.81) 



where as previously rij = rj — Tj and where we have added the extra (small) term 
in the denominator to prevent it becoming singular, this von Neumann-Richtmyer 
term becomes 



(nij)NR 



Pij 



where 



hvi 



f^ij 



(3.82) 



(3.83) 



By considering the bulk and shear viscosities of a generalised fluid it is possible to 
obtain a second form o: 



some time flLandshofl 



the artificial viscosity, a r id ind eed this has been known for 



1930 



Landau fc Lifshitz 



velocity divergence and uses the average sound speecc 
velocity component, giving the overall form 



19591 ). This form is linear in the 



Cs^ij as a second, characteristic 



(n.,)i 



asPH Cs,jj fJ'ij 
Pij 



(3.84) 



where /ijj is as before, and ctspH is a second constant of order unity. Note that the 
negative on the RHS arises from the requirements that the viscous force component 
must be non-negative (i.e. Uij > 0) and that it should be present only for convergent 
flows, where Vj^ ■ r^j < 0. In fact these criteria also hold for the von Neumann- 
Richtmyer form of the viscosity, and therefore in both cases the viscosity is set to 
zero in expanding flow conditions. 

These two forms have different and complementary numerical effects. At low 
Mach numbers (T W < 5) the linear form performs very well in shock tube tests 
f Monaghanl. 1985), whereas for stronger shocks it fails to prevent inter-particle pen- 



etration f Lattanzio et al 



19851 ). This is an unphysical phenomenon in which the 



two streams pass through each other at the shock front, leading to the possibility of 
two particles occupying the same position with differing velocities - a multi-valued 
velocity field. This possibility can be prevented by using the quadratic form of 



^Where as usual, the sound speed is defined as c^ = dP/dp. 
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von Neumann &: Richtmverl as it provides a stronger viscosity for high Mach num- 
ber, although conversely, on its own this decays too rapidly at low Mach numbers 
and fails to damp out the unphysical post-shock oscillations or "ringi ng" that oc- 



curs. The standard solution is therefore to use the sum of the two terms (JMonaghan 



19891 ) ■ resulting in a "standard" SPH viscous term of the form 



— aspH Cs,ijHij + ,5spH ^J^ij 
% = { T:, ^'^ ■ '^^' ^ ^' (3.85) 

otherwise. 

Various numerical tests have showed that in general the constant values aspH = 1, 
/3sPH = 2 and e = 0.01 in equation 13.831 give good results without significantly 
affecting non-shocked flows. However, throughout the simulations discussed in the 
later chapters of this thesis we have used values of asPH = 0.1 and /3sph = 0.2, which 
have been found to be adequate to accurately resolve (weak) shocks, while at the 
same time minimising the artific ial heating which wou ld have biased our simulation 
results - details can be found in 



Lodato fc Ricel (120041 ) 



This general form of the viscosity can then be incorporated into the momentum 
equation to give the following form; 






E ^^ (^ + ^ + n.*) ^^■^^■- (3.86) 



Given that the artificial viscosity term is also pairwise symmetric in i, j (since both 
Yij and Vjj are anti-symmetric in i and j) it is clear that this form of the equation 
of motion also conserves momentum exactly. Likewise it is clear that angular mo- 
mentum is conserved, and furthermore, by a similar argument to that presented in 
Section [3.3.31 it is possible to show that in order to preserve energy conservation, the 
energy equation must be modified to include an extra dissipative term such that 



^^ - ^ Y^ m.Vji ■ V^Wji + ^ niiUji-^ji ■ V^Wji. (3.87) 



,2 



dt p 

In this manner it is therefore possible to include a dissipative term such that shocks 
can be accurately captured, albeit broadened across a few smoothing lengths. Since 
mass, momentum and energy are still explicitly conserved across the shock the 
Rankine-Hugoniot equations are automatically satisfied at distances greater than a 
few smoothing lengths from the shock. Furthermore, since (in theory at least) the 
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artificial viscosity is zero away from shocks, all the other initial criteria are satisfied 
also. However, there are various improvements that can be implemented, and these 
will now be briefly discussed. 

3.4.2 More Advanced Viscosities 

The thorn in the side of all viscosity prescriptions is the requirement that in the ab- 
sence of shocks or other natural dissipative processes the artificial viscosity should 
reduce to zero, thereby requiring some means to discriminate between shocks and 
other flow features. Compounding the problem is the f act that careful consid eration 



of the artificial viscosity given above (see for instance iLodato fc Pricell2010[ ) shows 
that it provides both a bulk and a shear viscosity, while to resolve shocks only the 
bulk component is required. Any artificial viscosity in the form of equation 13.851 
therefore necessarily introduces an unrequired shear viscosity, which can be prob- 
lematic in situations where shear flows are important (such as discs), leading to 
spurious energy and angular momentum transport. Furthermore, since the shear 
force across any given particle varies with its smoothing length, it is clear that this 
shear component is resolution dependent. Generally speaking however thi s effect 



can be reduced by sensible choices for asPH and /3sph (JLodato fc Ricd . 1200J) - this 



further explains the low values of aspH and /3sph mentioned earlier. 



3.4.2.1 The Balsara Switch 



An at tempt to reduce the induced viscosity in shear flows was presented by [Balsara 



(119951 ). in which the standard artificial viscosity term IIjj is diminished by the factor 



fij = |/i + /jl/2, where 

■'' |V-Vi| + |V X v,|+0.0001cs,i//i' ^' ' 

The inclusion of the vorticity (the curl of the flow field) allows this form of the 
viscosity to perf orm better in shearing and obliquely shocked flows (see for instance 



Steinmetz 



19961 ). while remaining unaffected in the case of normal shocks. In a 



similar manner to the "standard" artificial viscosity term, this form also includes a 
small term O.OOOlCg j//i to prevent the viscosity from becoming singular. 
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3.4.2.2 The Morris & Monaghan Switch 

Although the Balsara switch represents a considerable improvement over the stan- 
dard form of artificial viscosity, problems still arise in the case of shocks in shearing 
flows, such as those found in accretion discs. For this reason, iMorris fc Monaghan 
(Il997l ) introduced the idea of a time-variant viscosity such that IIjj remains un- 
changed from the standard form, but where asPH = a(^), and where /3sph = 2q;sph- 
The value of a is then evolved for each particle according to the following equation; 



da 

d^ 



a 



Or 



+ S^. 



T 



(3.89) 



Here amin ~ 0.1 is some minimum value, justifled by the requirement that some 
level of artiflcial viscosity is required to maintain particle ordeio, r ~ 0.1 — 0.2 h/cg 
is a decay timescale (chosen so that the viscosity decays away over a few smoothing 
lengths) and S^ = max(— V • v, 0) is a source term, activated whenever the flow 
becomes convergent. Although this form of the source term is still non-zero for pure 
shear flows, this is counter-balanced to some extent by the decay term, and ha s been 
found to work well in many tests of the artiflcial viscosity (JDolag et al.l . l2005l ). Fur- 
ther variations on this theme have been effected, including incorporating the Balsara 
switch into the Uij term, and capping the maximum value t o which a can rise b y 
using a source term of the form S^ = max((amax — a)V-v, 0) (JRosswog et al.l . l2000l ). 
For a good general overv iew of the relative mer i ts of a va r iety o f artiflcial viscosity 
methods, see for instance 



Lombardi et al. 



fll999h : lRosswoelf!2009f ): 



CuUen fc Dehnen 



(|2010| . in prep. 



3.4.3 A Note on Entropy 

All of the above methods have essentially been aiming to capture the same phe- 
nomenon, namely the increase in entropy found across a shock front, while simul- 
taneously ensuring isentropic flow elsewhere. Furthermore all share the common 
feature that flow evolution pro ceeds via integration o f the energy equation. How- 
ever, an approach espoused by ISpringel &: HernquistI (120021 ) is to consider evolving 
the entropy directly, thereby ensuring that the entropy can only increase. 

In this manner, we recall that in terms of density p and speciflc entropy s, the 

0.01 should 



^No te that only very low lev els of viscosity are required for this purpose - amin 
suffice (jCuUen fc Dehnenl . l2010l in prep.). 
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equation of state is given by 



P^ = t^i{si)pl (3.90) 



for some entropic function k{s). Similarly, the internal energy Ui may be obtained 
from p and s via 

u. = ^p]-\ (3.91) 

7-1 

In the case of isentropic flow, we have k{s) = const, and thus by definition 

CI rbT , , 

However, in the case where artificial viscosity is included, the time derivative of the 
entropic function becomes 

dn, 1 7 -1 J2m^U,J^r,J ■ V.W^j. (3.93) 



dt 2pJ 



By noting that 



V,iy,, = \ViWij\rij, (3.94) 



and also that IIjj is only non-zero for Vjj ■ Vij < 0, it is clear that the term on 
the RHS of equation 13.931 is strictly non-negative, and thus that entropy can only 
increase throughout the flow. Using this method of evolving the flow properties it 
is therefore possible to explicitly ensure that the entropy of any particle increases 
monotonically with time. 

3.5 Variable Smoothing Lengths 

Up to now, it has been assumed that the smoothing length h is held constant with 
time, and is moreover equal for all particles. In regions where the density (and thus 
the number of neighbours) is roughly constant, this maintains a constant (small) 
sampling error within the SPH smoothing kernel. This requirement of constant 
smoothing length is quite restrictive however, as it prevents the code adapting ef- 



fectiv ely to regions of higher or lower than average density (jSteinmetz fc Mueller 



19931 ). By allowing the smoothing length to vary both temporally and spatially, sam- 
pling errors can be minimised across regions of varying density, as either the number 
of neighbours or the mass within a smoothing kernel (and thus the resolution) may 
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be maintained. There are various simple ways of all owing: v ariab le effective smooth- 
ing lengths that have been introduced, for instance iBena f ll990l ) suggested using a 
symmetrised smoothing length hij = {hi + hj)/2, such that the kernel becomes 



W., 



w 



hi + hj 



U' 



(3.95) 



An alternative method has been suggested by iHernquist &: Katzl (119891 ). in which 
the average kernel value is used rather than the average smoothing length, such that 



W, 



W{rij,hi) + W{rij,hj 



u 



(3.96) 



With variable smoothing lengths it then becomes necessary to determine the value 
of h for each particle. A standard method of doing this is to link the smoothing 
length to the local density, such that 



Pih! 



const. 



(3.97) 



Since this constant clearly has units of mass, it is frequently linked to the particle 
mass, giving the following prediction for the particle smoothing length; 



hi = T] [ — 
Pi 



1/3 



(3.98) 



wher e the coupling constant is generally in the range 1.2 < r] < 1.5 ( jRosswogi . 
20091 ). By construction this method maintains a constant mass within the smoothing 
kernel. As each of the above formalisms remains pairwise symmetric, momentum 
remains fully conserved, and increased spatial resolution is achieved at relatively 
low computational cost. The latter method (using the averaged kernel value as in 
equation l3.96l) has additional advantages in it is less problemat ic across shocks, and it 
coup les better with tree methods for calculating self-gravity (jSteinmetz &: Mueller! . 
19931 ). Nonetheless, in both cases errors appear in either the entropy or energy 
equation, such that either 



dE dn(s) 

— — or — ; — 

dt dt 



dWdh 
'dh'dt 



7^0, 



(3.99) 



(JHernquistl . Il993l ). and the relevant quantity is therefore not explicitly conserved. 
It is however possible to construct SPH estimates that self-consistently account 
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for the variation in smoothing length, and therefore ensure exact energy conserva- 
tion. In this case, the estimator for density equation 13.391 becomes 

p. = ^m,iy(r,,,/i,), (3.100) 

i 

noting that the smoothing length used in the kernel is that associated with particle 
j only, and thereby remains constant throughout the summation. By taking the 
(Lagrangian) time derivative, we obtain 

^^^™.(v„.V<W„(.,,.fif). (3.101, 



noting the extra terms compared to equation 13.411 and where now we set Wji 
WijCji, hj). Noting that 

(3.102) 



dhj dhj dpj 



dt dpj dt 
and using equation 13.981 we see that 

dhi hi 

-^ = --^. 3.103 

Substituting this into equation 13. 1011 and gathering like terms, we find that the time 
derivative of the density becomes 

J A 



where 



hi v-^ dWij(hi) 

and where dWji/dhj is known from the choice of kernel. Although it can be cal- 
culated directly from the kernel, in the case of the cubic spline kernel given in 
equation 13.311 it is generally evaluated by noting that 

dW 3 

= -xVW - -W, (3.107) 

dhj h 
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where W and VW are given by equations 13.311 and 13.331 respectively. 

Similarly, there is a correction factor to the momentum equation to allow for 
the spatial variation in smoothing lengths. Recall from equation 13.521 that in order 
to calculate the spatial variation of the Lagrangian, we need to know the spatial 
derivative of the density. Allowing now for variable smoothing lengths and using 
equation I3.57| we therefore find that 



E™.(v,w,(.,)fe-M.«l^t)- '"-' 



By gathering like terms, we find that the correction factor for the spatial derivative 
of the density is same as that for the temporal one, namely that 

^ = 7^ 5Z"^fcV,l^,fc(/i,) [6,i - 6,k] , (3.109) 



■' k 

with the factor ^2^ defined as before in equation 13.1061 

Following the same derivation as in Section l3.3.2| it is then easy to show that the 
acceleration due to hydrodynamic forces with spatially varying smoothing lengths 
is given by 

dv / P P \ 

Finally, from equation I3.73| we see that the evolution of the internal energy in the 
presence of variable smoothing lengths becomes 

^=(52E"^^^^-^-^^^-^(M- (3.111) 

J' J i 

By an analogous process to that described in Section l3.3.3[ it is possible to show that 
this equation for the evolution of the internal energy is also explicitly conservative of 
the total energy of the system, E. The three equations 13. 100| 13.1 101 and 13. 11 II along 
with the relationship between the density and the smoothing length equation 13.981 
therefore form a fully consistent, fully conservative SPH formalism with spatially 
varying smoothing lengths. 

A problem exists however, in that in order to obtain the density, one needs to 
know the smoothing length (equation 13. lOOp and to obtain the smoothing length one 
needs to know the density (equation 13. 98p . In order to resolve this, this pair of equa- 
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tions can be solved iteratively by the Newton- Raphson method (jPrice fc Monaghan 



20071 ). By rewriting equation 13 .981 we can combine these two equations to reduce 



the problem to that of finding the root hj of the equation ({hj) = 0, where 

Cihj)=m,(^^^ -^m,iy(r,„/i,). (3.112) 

Here the first term represents the density obtained from assuming a fixed mass 
within the smoothing kernel, while the second term is the standard SPH estimate 
for the density. From some initial estimate of the root hj, the Newton- Raphson 
method gives a better estimate as being 



CM 



lijjnew — lij /-iCL \i (^d.iidj 



where the prime denotes differentiation with respect to h. By using equation 13.1061 
we see that 

ah,) = -^, (3.114) 

and thus the updated value /ij,new is given by 

This may be repeated until |/ij,ncw — hj\/hj < e for some small value of e, frequently 
set to 10^'^. Then in turn, a self consistent value of the density is then obtained from 
equation 13.981 As there is generally relatively little change in hj and pj between 
timesteps, the estimator for hj is taken as the value from the pr evious timestep, and 



conve rgence usually occurs within a small number of iterations (jPrice fc Monaghan 



20071 ). In pathological cases where the Newton- Raphson method does not converge, 
other, universally convergent but slower methods such as the bisection method may 
be used instead. 

Although the inclusion of variable smoothing lengths through this method does 
inevitably increase the computational cost of the code, this is relatively small, and 
the conservation properties are recovered to within machine (and integrator) tol- 
erance. Other tricks, such as predicting the change in the smoothing length using 
equation 13.1021 can redu ce the computational cost still further (see for instance. 



Price fc Monaghanll2007f ) 
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3.6 Including Gravity 

As many astrophysical situations are driven at some level by gravitational forces, it 
is important to be able to include this consistently within the SPH framework, and 
in such a manner that the inherent conservation properties of the algorithm are not 
compromised. While much work has been put into N-body simulations of discrete 
particles, within the SPH formalism we are aiming to model the gravitational force 
over a continuum, and thus it should be smoothed (or softened in SPH parlance) in 
a similar manner to that in which the discrete particle mass is smoothed into the 
density of a fluid continuum. In this section we therefore consider how this can be 
done in a consistent manner, and one in which as before momentum and energy are 
explicitly conserved. 

3.6.1 Gravity in the Lagrangian 

In an extension to the Lagrangian for the hydrodynamic equations of motion, it is 
possible to incorporate the effects of gravity by considering a Lagrangian of the form 

£(r,v) = ^miQvi-v,-M,(ri)j - *, (3.116) 

where \1/ is an as yet undefined measure of the total gravitational potential energy of 
the system. By comparison with equation 13. 46l this is clearly just the hydrodynamic 
Lagrangian with an additional term 

>Cgrav = -^ (3.117) 

which describes the effects of gravity. 

Now, as with the density, at position Fj we can obtain the local gravitational 
potential $i, via a sum over all particles such that 

$i = G'J]m,0(r,-r,-,ei), (3.118) 

j 

where 0(rj — r^) = (pi^Vij) is known as the (gravitational) softening kernel, G is the 
universal gravitational constant and where Ei is the softening length associated with 
particle i. The softening kernel at this stage is fairly general, but it must have the 
following properties: 



30 



Smoothed Particle Hydrodynamics 3.6. Including Gravity 

• 0(r, h) < for all r, h, as the local potential $ must be strictly negative 
definite; 

• V0(O, h) = 0, such that the gravitational force exerted by any particle on itself 
is zero; 

• lim </'(r, /i) = — , i.e. the softening should reduce to zero at large inter- 

r/h->-oo r 

particle distances, and the Newtonian potential should be recovered. 

Generally speaking, and throughout this thesis, it is assumed that the softening 
length is exactly equal to the smoothing length for all particles, i.e. Si = hi for all 
i. In a similar manner it is generally taken that the force should only be softened 
when r < 2h, so that force softening and density smoothing occur over exactly the 
same region. 

Noting that the gravitational potential energy is just the mass times the gravi- 
tational potential, since the latter is defined over pairs of particles, by definition the 
total gravitational potential energy of the system is given by the sum over all pairs 
of particles, such that 

^ = G^mi^mj(j)ij{hi) (3.119) 

i j<i 

= -^^^^i^j^i:i(^i) (3.120) 

i J 

Note that equation 13.1191 sums over all pairs of particles, including the so-called self- 
interaction terms where i = j, and thus explains the factor of a half in equation l3.120l 
From this definition we therefore find that 

^ = ^$^rn,<D„ (3.121) 

i 

and thus that the full Lagrangian in the presence of gravity becomes 

£(r, v) = ^ m, Qv, ■ v, - M,(r,) - i$, j . (3.122) 

By considering only the gravitational term in the Lagrangian, we can as before 
use the Euler-Lagrange equations 13.471 to obtain the acceleration due to gravity. 
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which becomes 

™,^ - ^. (3.123) 

Using equations 13.1171 and 13.1201 we therefore find that the spatial derivative of the 
gravitational Lagrangian becomes 



-d^ = -^1^1^^^^^^^ (3.124) 

■^ i k •> 

i k ^ * J / 

Here we see that in the case of fixed smoothing (and therefore softening) lengths, 
dhi/drk = 0, and thus we only require the first term to determine the effects of 
gravity. The second term is therefore a correction term to allow for spatial variation 
in h. 

As before, using the method of equations 13.541 to 13.571 the spatial gradient be- 
comes 



L^-'-'grav 



dVj 



G . 

— ^ ^ minik VjCpikihi) [6ij - 6kj] , (3.126) 



2 

i k 



G - Q 

" 9""^J X] i^k^j(t)jk{hj) + IT X] mimjVj(f)ij{hi). (3.127) 



Now by changing the summation index of the first term to i, and noting again that 
the kernel is antisymmetric we obtain 



^^-grav 



dr, 



G 

-jTT^j Yl ^^ i^Ajii^j) + ^Ajiihi)) , (3.128) 



which therefore encapsulates the effects of gravity in the case of constant smoothing 
lengths. 

If we now consider the second term in equation 13.1251 and self-consistently correct 
for spatial variation in the smoothing length, we find that 



djC 



dVj 



i k ■' 
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By substituting equation 13.1091 into the above we find that 






grav 



Or, 



= - f E E "<'"' ^^SS- E '"'^^W'"(AOfe - *il (3.130) 

i k I 1^1 I ^ 

i k 

Now by changing the summation index of the second sum in the first term of equa- 
tion [3TT3T] to i, defining a new quantity ^p such that 

«.-|;E™,^. (-32) 

and using the antisymmetry property of the gradient of the smoothing kernel, we 
see that the correction term reduces to 



dr. 



dvj 



~m,J2^, (^V,iy,.(/i,) + ^V,W,,{hA . (3.133) 



Finally, using equation 13 . 1 231 and by incorporating the effects of gravity into the 
equations of motion for a hydrodynamic flow with artificial viscosity (while self- 
consistently allowing for variable smoothing lengths) we find that the full equations 
of motion become 

^ 5^m, (V,-0,,(/i,) + V,0,,(/i,)) (3.134) 



dt 

~ 2 



T E^^ (^V,iy,.(/i,) + ^V,W,.{h,)^ 



with fij and ^i defined as per equations 13. 1061 and 13. 132] respectively I. As in the case 
for the pure hydrodynamic flow, the use of a Lagrangian in deriving these equations 
ensures the explicit conservation of both linear and angular momentum, which is 
also clear from the pairwise symmetry present in all terms in the above equation. 



^Note that for consistency, the artificial viscosity term Hji uses the average value of the smooth- 
ing lengths hji ~ {hj + hi)/2 in its definition of fiji (equation 13. 83p . 
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3.6.2 Evolution of the Gravitational Potential 

Clearly as particles move about within a gravitational potential, their potential 
energy (given in SPH terms by rrij^j) will also vary. Although the potential (and 
thus the potential energy) is obtained at any point by the sum over particles using 
equation I3.118| the time evolution of the potential energy is required to maintain 
energy conservation. Hence in a similar manner to Section 13.3.31 we must consider 
the total energy of the system, which including the gravitational potential energy 
becomes 

^ = 5Z ^J- [l^J ■ ^1 + % + l^j) ■ (3-135) 

j 

As before, to ensure energy conservation we require that the time derivative of the 

total energy is zero, i.e. that 

E/ dv, dui 1 d$A 

J 

By considering equation 13.1181 we see that 



d$j _ G 

"df ~^ 



E- (V,^.(y ■ ^ + V.^,(/„ . ^ . ^^^ . (3.137) 



Recalling the definition of ^j from equation I3.132| and using equation I3.1U4I for the 
definition of dpj/dt with variable smoothing lengths, we obtain 



d<l>j _ G 
"df ~ ^ 



J2 m,^r^^ " (^jM^j) + ^,^,W,,{h,)\ , (3.138) 



From Sections 13.3.31 and 13.51 we know that in the energy balance (equation 13. 135p . 
over all particles the hydrodynamic terms in the equations of motion (equation 13. 11 Op 
exactly counteract the temporal rate of change of the internal energy. 



dE'hydro _ ■^;^ ^ ( ^^ dvj 

hydro 



^2"^^ i^i 



dt ^^ M ■' dt 



^1=0, (3.139) 



and thus pure hydrodynamic flows are exactly conservative of energy. With the 
inclusion of gravity we therefore only need to show that over the whole system the 
gravitational terms in the equations of motion balance the time derivative of the 
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gravitational potential, i.e. that 

ld$ 



(i^^grav _ >r^ / dvj 



Y. ^i [ ^i 



dt ^^ M ■■ dt 



J 



grav 2 dt 



(3.140) 



in order to maintain exact conservation of energy in self-gravitating systems. 
From equations 13.1341 and 13.1381 this gravitational energy balance becomes 



dE 



dt 



G ( 



^V,W,,(/i,) + ^V,iy,.(/i.) ) (3.141) 



G [ £ 



Cancelling like terms, this reduces to 



dE 

dt 



+ ^v,-v,iy,.(/i.) + |-v,iy,.(/.,)) 

^ 3 / 



(3.142) 



and finally, noticing that the gradients of both the smoothing and the softening 
kernels are antisymmetric under a reversal of the summation indices i and j, we 
obtain the desired result that 

^^ = 0. (3.143) 

Therefore, we see that gravity can be included into SPH in such a manner that the 
algorithm remains explicitly conservative of energy. 

3.6.3 Gravitational Potentials and the Softening Kernel 

Finally for this section, we need to consider the form of the gravitational softening 
kernel 0, and its relation to the smoothing kernel W . Recall that Poisson's equation 
links the gravitational potential $(r) to the density p(r) at position r, such that 

V2$(r) = 47rGp(r). (3.144) 
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Given that we implicitly assume each particle to be spherically symmetric, by using 
spherical polar co-ordinates and substituting equations 13.391 and 13.1181 into equa- 
tion 13.1441 we find that (for a generalised radial co-ordinate r) 



W{r, h) 



d 



,d(l){r,h) 



(3.145) 



Airr"^ dr \ dr 

where we have neglected the spatial variation of /£. 

We can now integrate this, to link the derivative of the softening kernel d(j)/dr 
(also known as the force kernel) to the smoothing kernel, such that 



dr 



4:71 



/r 
r'^W{r')dr' + 



9l 

^2 ■ 



(3.146) 



with the integration constant Ci subject to the condition that for r > 2h we recover 
the standard Newtonian inverse square law, which using our definitions becomes 
d(f)/dr = 1/r^. In a similar manner we can integrate this a step further (by parts), 
to give the full softening kernel, such that 



An 



1 



r'W(r')dr'+ / r'W{r')dr' 






(3.147) 



where the second integration constant allows the correct asymptotic behaviour (i.e. 
— )■ as r — 7- oo) to be established. 

With this in mind, for the cubic spline kernel defined in equation 13.311 the force 
kernel d(j)/dr becomes 



(90(r, h) 
dr 






6 3 1 4 

^ 



63 1 4 

5 6 



15a; 






<a: < 1, 

1 < X < 2, 

x> 2, 



(3.148) 



V. T 



where x = r/h and where the integration constants have been absorbed to ensure 
piecewise continuity. Finally, we therefore find the full softening kernel consistent 



^This is because the smoothing length essentiahy acts as a normahsing constant in both the 
smoothing and the softening kernels, and for any given particle is held constant within Poisson's 
equation. Thus its spatial variation is immaterial here. 
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with the cubic sphne smoothing kernel to be 



0(r, h) 




3 3 1 5 7 
10 10 5 



X 



30 



-X 



15a; 



<x< 1, 

1 < X < 2, 

x>2. 



(3.149) 



Using this definition of the softening kernel along with the cubic spline smoothing 
kernel equation 13.311 the equations of motion equation 13. 1341 and equation 13. 1381 for 
the evolution of the gravitational potential therefore allows gravity to be included in 
a manner that it is fully conservative, and is such that Poisson's equation is satisfied 
throughout. 

3.7 Finding the Nearest Neighbours 

Various methods exist for finding the nearest neighbours (i.e. those particles within 
the smoothing kernel of any given particle), with the simplest being a direct search 
over all particles. This becomes very expensive in the limit of large numbers of 
particles A^ however, as the computational cost scales as O^N"^). Other methods 
such as using an overlaid grid or a linked list of partic l e positions have bee n used 



f Hocknev fc Eastwood 



1981 



Monaghan 



1985 



Murrav 



1996 



Deeganl . l2009h . One 



of the more efficient methods however is to use a hierarchical tree structure, an 
approach that grew out the requirements of N-body codes to distinguish distant 
particles (where the gravitational forces could be evaluated via multipole expan- 
sions) from local particles (where direct N-body calculation of the forces was still 



lese in genera l reduce the cost of nei 



required) . T 

C(iVlogA^) fiBarnes fc Hut 



1986 



Hernquistl . Il987 : 



ibour-finding from OiN -^) to 



Hernguist fc Katzl . Il989l ). al 



Dehnen 



2002 



2000|). 



though reductions to 0{N) have been achieved (|] 

Trees are essentially data structures which decompose the computational domain 
into a series of discrete volumes, the sum over which contains all the particles. The 
smallest volumes generally contain only a single particle, but this i s not strictly 
necessary, and indeed may inhibit the efficiency of the code (JDehnenl . l2009l . private 



communication). By construction, particles which are near each other in space are 
near each other within the tree structure, and thus by looping over a relatively small 
part of the tree, the nearest neighbours may be found efficiently. There are various 
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different algorithms that perform this decomposition described in the hterature, none 
of which are trivial, so I shall not at t empt t o go into any d e pth here. For further 
details, see for instance 



19930; 



Barnes fc Hut 



Press et al. 



feoOTri: behnenl ( 120021 ): 



Steinmetz Sz Mueller 



(119861 ): JBentlevi (119751 ) and references therein. 
A distinct advantage of using trees for neighbour finding within an SPH code is 
that they couple readily with pre-existing methods for evaluating the gravitational 
force between large numbers of particles. Rather than a direct summation (which 
is of 0(A^^)) over all particles to find the gravitational force at a specific location, 
particles at large distances can be effectively treated as a single body, and multipole 
expansions used to approximate the force. This approach has found success in 
various A^ -body codes as a means of reduc i ng the computation time to 0( N log A^) 



or lower ( iBarnes fc Hut 



1986 



Hernquistl . 119871 : iHernquist fc Katz 



1989). Use of 



a tree algorithm therefore allows the process of neighbour finding to be coupled to 
that of finding the gravitational forces acting on a particle, with an attendant saving 
of computational expense. 

3.8 Integration and Timestepping 

So far we have obtained equations to evolve the density, the three components of 
velocity under the influence of pressure, gravitational and (artificial) viscous forces, 
the internal energy and the gravitational potential. Finally therefore it is time 
to consider how these equations are actually evolved, and to discuss the issues of 
temporal integration and time-stepping. 

Generally speaking there are two principal methods used to perform the time 
evolution, and indeed the code I have used throughout includes the option to use 
either. They are the so-called Leapfrog integrator (also known as the kick-drift or 
Stormer-Verlet integrator) and the Runge-Kutta-Fehlberg method, and I shall now 
briefly consider both of these. 

3.8.1 The Leapfrog Integrator 

The leapfrog integrator is a second-order integrator, so-called because the position 
and the velocity are advanced half a timestep out of phase, with each update of 
position or velocity using the value of the velocity or position evaluated at the 
previous half timestep. In this manner, the positions and velocities "leap-frog" over 
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each other at every half timestep, giving rise to the name. The leapfrog method is 
widely used in A'^-body codes, since in the case where the acceleration is independent 
of the velocity, i.e. where a = a(r) only, it is particularly simple to implement. In 
its "pure" form it is a time-reversible, symplectic integrator, which by definition 
is explic i tly co nservative of both energy and angular momentum (see for instance 



Springell ( l2005l ) and references therein). 

In essence then, if the position, velocity and acceleration at time tj are given 
by Fj, Vj and a^ respectively, with a timestep 5t the standard form of the leapfrog 
integrator gives the positions and velocities as 

' (3.150) 

Vi+i/2 = Vi_i/2 + aj bt. 

Here it is clear that the positions and velocities are evaluated at half timesteps with 
respect to each other, and "leap-frog" over each other as they are evolved. In this 
form it is also clear that the integrator should be perfectly time-reversible. 

A form that is often more readily applied is the equivalent definition at integer 
timesteps, which becomes 

. ( St \ 

Tj+l = Fi + c)t Vi + — aj , 



V ^ / (3.151) 

ot 
Vi+1 = Vi + — (ai + a^+i) . 

Although it is now less obvious, these equations are still fully time reversible. Note 
further that the form of the increments on the RHS of each of the above equations 
is equivalent to an estimate of the relevant quantity at the following half timestep, 

noting in particular that 

5t 

Vj + — ai = Vi+i/2, (3.152) 

to first order. 

Note however, that in both cases problems arise if the acceleration depends on 
the velocity, since from equation 13.1511 we see that to calculate Vj+i we already 
need to know the acceleration a^+i, and the scheme becomes implicit. Since in 
SPH simulations both the pressure force and the artificial viscous force depend on 
the local velocity, it is clear that modifications are required before this integrator 
may be used. The standard way this cor rection is implemented (see for instance 



Springel et al.ll200ll : IWetzstein et al.ll2009l ) is as follows: 
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Firstly, predict the positions at time tj+1/2 in a manner analogous to equa- 
tion 13.1521 via 

ri+1/2 = ri + — Vi. (3.153) 



Secondly, use equation 13. 1521 to obtain the velocity at time ti+1/2, and extrapo- 
late other values (such as density, internal energy and gravitational potential) 
at the half timestep also. Hence calculate the acceleration at the half timestep, 

aj+l/2- 



Calculate the velocity at time tj+i using 



'i+\ 



6t Sij 



i+l/2- 



(3.154) 



Now update the positions to timestep tj+i using 

6t , 



' j+i 



Vi 



'i+i 



)• 



(3.155) 



The process is now repeated as required. 

Although the strict symmetry between the integration of positions and velocities 
has been lost by the inclusion of these predictor steps, this method still remains 
time-reversible. Furthermore, it is also now possible to include adaptive timestep- 
ping, which would have been problematic before pr ecisely because of the symme- 
try between the integrations (jWetzstein et al.l . l2009l ). Generally spe aking however, 



mairi t aining time-reversib ility with adaptive timestepping is difficult (jSpringel et al. 



2001 



Quinn et al.l . 119971 ). though not impossible. 



3.8.2 The Runge-Kutta-Fehlberg Integrator 

Runge-Kutta methods for integrating systems of differential equations are well 
known, tried and trusted methods, which use multiple estimates of the derivative 
across a given timestep to arrive at accurate, generally high order estimates for the 
evolved quantities. Most common is the fourth order Runge-Kutta method, often 



simply abbre viated to R K4, which has 



flKuttal . Il90lh . Moreover, iFehlberd ( 11968 



j een k nown and used for over a century 



19691 ) obtained a modified Runge-Kutta 



integrator (now known as a Runge-Kutta-Fehlberg integrator) which embedded a 
order n + 1 method within an order n method. This allows the two methods to be 
compared to give an error estimate, and thus for the error to be controlled to some 
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given tolerance. The most common of these methods embeds a fifth order estimate 
within a fourth order scheme, and is therefore known as an RK45 integrator. 

However, compared to the leapfrog integrator, which requires only one evaluation 
of the acceleration per timestep, the RK4 scheme requires four, and the RK45 
method requires six. Therefore these methods, although correct to much higher 
order than the leapfrog, add significantly to the computation time required. (Note 
also that they are not necessarily more accurate either, as they are not explicitly 
conservative in the way that the leapf r og ra ethod is. See for instance ISpringel 
(120051 ) ; IWetzstein et al.l (120091 ) ; iRosswod (120091 ) for a comparison of these integrators 
as applied to a simple Keplerian orbit evolved over many dynamical times.) The 
solution is to go to a lower order RKF method, where the implicit error control is 
still present but the number of derivati ve evaluations is r educed. A common choice 



for many SPH codes including VINE (IWetzstein et al. 



2009 



Nelson et al. 



20091) 



and the one us ed in the code I have used, is the RK12 integrator developed by 
FehlbergI (11969! ) which proceeds as follows. 

For a given variable z, the evolution from Zj at time tj to z^+i at time tj+i = ti + 6t 
is given by 



'H+l 



1 



255, 



h St, 



(3.156) 



256 " 256 

where the values of ko and ki are provided by evaluating z at various points, such 
that 



ko = z(ti, Zj), 

6t 6t, , 

ki = z{ti + — , Zj + —ko] 



(3.157) 



and where the dot as usual denotes differentiation with respect to time. Expan- 
sion via Taylor series shows that this is accurate to first order, with the choice of 
coefficients in equation 13.1561 producing a leading order truncation error Ttrunc such 
that 

nmnc = -777:<^i^z- (3.158) 



512 



Using the values for ko and ki defined above, we can compute a further estimate 
z*_,_j^ for z at time tj+i, such that 



6t 



^i+l ~ '^i 



1 , 255, 1 , 

"^0 + 7:r;:ki + -^-^kr 



2 V512 " 256 512 



(3.159) 
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with the additional value /c2 defined such that 

/^ 1 , 255, \ , 
k2 = z{ti + dt,Zi+[-—ko + —-h]6t), 

V256 256 J (3.160) 

= z(tj+i,Zi+i). 

Again, by considering Taylor series expansions, this value z*_^_^ can be shown to be 
a second order estimate. One of the more appealing tricks of this method is that 
here /c2 is simply /cq evaluated for the next timestep, and thus per timestep, only 
two derivative evaluations are required. 

We now therefore have both a first and a second order estimate for the value of 
z at time tj+i, with a known truncation error for the first order method. This can 
therefore be used for error control, to ensure that the timestep used is appropriate 
(see for instance, 



Press et al. 



20071 ). However, in order for this error control to be 
valid, the first order scheme must be used for the evolution. To mitigate this, by con- 
struction this first order scheme has very small second order errors (equation [STTSH]) , 
and so is effectively a quasi-second order integrator. 

3.8.3 Timestepping Criteria 

For either integrator, it is crucial that the timestep size is chosen correctly, both 
to ensure the accuracy of the evolution and to ensure numerical stability. In this 
section I shall briefiy discuss the principal timestepping criteria in general use, and 
one specific to the code I have used. 

3.8.3.1 CFL Criterion 

By far the most general timestep criterion for gas-dynamical systems is the so-called 
Courant-Friedrichs-Lewy or CFL condition, given in it simplest form by 

ox 
StcFL < — , (3.161) 

c 



where Sx is a characteristic length scale, and c is a characteristic speed (jAnderson 



19951 ). For SPH simulations, these are both well defined; the smoothing length h 



provides the characteristic length, and sound speed Cg gives the characteristic speed 
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of the medium. The CFL condition for particle i then becomes 

OtcFL < —■ 



(3.162) 



This has a ready physical interpretation in that it prevents spatial information trans 
fer through the code at a rate greater than the local sound spee d. In the presence 
of artificial viscosity this requires a slight alteration, and as such 
recommend using the following criterion; 



Bate et al 



fll995[ ) 



St 



0.3/1 



CFL 



c, + /i| V ■ v| + 1.2(asPHC. + /3sPH/i| V • v|) ' 



(3.163) 



where the factors of 0.3 in the numerator and 1.2 in the denominator are empirically 
determined. The aspH and /3sph terms are those used to determine the strength of 
the artificial viscosity (see Section [33]), and it should be noted that the final term 
in the denominator is only included in the case where | V ■ v| < 0. The extra /i| V ■ v| 
term in the denominator accounts for the expansion or contraction of the flow, and 
thus explicitly al lows for compressibility effects . There are variat ions on this theme 



(see for instance lDeeganll2009l : iMonaghan 



1992 



Monaghanlll989l ) but the definition 



given above is the one present in the code I have used. 

3.8.3.2 Force Condition 

A further commonly used timestep condition is that based on the acceleration of 
the particle, known as the force condition. This is simple in form, and is given by 



6U 




(3.164) 



where as before a is the particle acceleration, and /p < 1 is a tuning constant. 
Values for fp vary from code to c ode but are generally in the range 0.25 - 0.5 



(IWetzstein et al. 



2009 



Bate et al 



19951 : iMonagharu . Il989l ). The code I have used 



employs /p = 0.3. 

3.8.3.3 Integrator Limits 

Dependent on the choice of integrator, other timestep criteria may be required. In 
particular, if using the RKF method the timestep criterion associated with the error 
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correction must be incorporated. Using the method outhned above, this corresponds 
to a timestep of 

StRKF = Sto\d\ -. 7, (3.165) 

y \'ZRK2 — ZrkiI 

where e is the desired error tolerance (usually of the order of 10~^ - 10~^) and the 
'Zrki,'Zrk2 are the predictions for any quantity z from the first and second order 
methods within the integrator respectively. The 5toid term is simply the increment 
used for the previous step. 

3.8.3.4 Generalised Timestep Criteria 

A general class of additional timestep criteria may be obtained by dimensional anal- 
ysis, in that for any time- varying quantity z we may define a characteristic timescale 
on which it varies as 

r. = -, (3.166) 

A/ 

where as usual z is the time derivative of z. To ensure that this timescale is properly 
resolved, we can therefore define a timestep condition such that 

5t. = /.|, (3.167) 

where /^ < 1 is a tuning factor. Although seldom required in general, a timestep 
criterion of this form was implemented into the code when looking at the effects of 
strongly varying cooling times in Chapter 5, and will be discussed in more detail 
there. 

3.8.4 Setting the Timestep 

There are therefore a variety of possible timestep choices, and thus to ensure that 
they are all satisfied, the timestep for each particle used is the minimum of all 
possibilities, i.e. 

^ti = min{6tcFL,i, Stp^t, (5tRKF,i, Stz). (3.168) 

Where there are only relatively small changes in the characteristic timescale, a stan- 
dard choice is to use a global timestep ^tgiob, which is set by the minimum of the 
timesteps Sti for the individual particles, such that 

Stgioh = min((5ti). (3.169) 

i 
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This has the advantage that all particles are evolved in lockstep, and thus there is 
no 'information lag' due to particles being on separate timesteps. 

On the other hand, individual particle timestepping has the advantage of be- 
ing much faster and thus more computationally efficient wherever there are large 
ranges in the timescales of the probl em being investigated. It can however intro- 



duce instabilities into the integrator f lWetzstein et al.l . |2009| ). and can also lead to 



the phenomenon of low density particles on long timesteps drifting into regions of 
high den sity evo l ving on much shorted timesteps, leading to spurious entropy gen- 



eration (jPearcd . l20ld . private communication). Th is latter effect is p articularly 



noticeable in tests of Sedov blasts (see for instance iTasker et al.ll2008l ). in which 
small entropy-driven bubbles lead to granularity in the post-shock region. This is 
a relatively uncommon phenomenon however, and occurs principally in the case of 
strongly shocked systems. 

Integrator stability may be maximised (particularly for the leapfrog scheme) by 
using timesteps that are integer multiples of each other, and generally speaking, for 
a maximum timestep T, sub-timesteps will be given by 2~"r. The particle timesteps 
are therefore rounded down to the nearest relevant power of two in this case. This 
is the case in the code I have used for all the simulations presented in this thesis, 
which uses individual particle timesteps and is only weakly shocked throughout. 

3.9 Sumraary 

In this chapter I have derived the SPH algorithm from first principles, and then 
built it up in a series of steps to solve for pure hydro dynamical isentropic fiows, dis- 
sipational fiows, and finally dissipational fiows under the infiuence of gravitational 
forces. Additionally I have shown that it is possible to self-consistently allow for spa- 
tially variable smoothing lengths, which allows the algorithm to be highly adaptive 
with the local fiuid density, but to maintain exact conservation of mass, linear and 
angular momentum and energy, to within the integrator tolerance. In the case of 
isentropic fiows, entropy is also conserved by construction. Furthermore I have also 
briefiy detailed various methods of finding the nearest neighbours, and two means 
of evolving the fiuid fiow forward in time. 

Since the problems I shall be considering in later chapters require only that 
dissipational fiow in the presence of gravity to be modelled, this is all that I 
have covered here. However this is by no means the limit for the SPH for- 
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malism. Much effort has been put int o including; additiona l phys ic s such as ra- 
diativ e transfer flNavakshin et al.l (120091): iPetkoya fc Springell fl2009f): iForgan et al 



(see for instance 



20091) :lBisbas et all (l2009f ): lGritschneder et all ( l2009h : lPawlik fc Schavd (l2008r ) and 
Altav et al.l (l2008f) to name b u t few of the recent e f forts) and magnetic field s /MHD 

Dolag fc Stasvszvnl J2OO9I): Irdsswoet fc Pric3 J2OO7I ): 



Price 



201 



i. 



2004 ) ■ and this will no 



Price fc MonaghanI (J2005l . |2004| ) and iPrice fc Monaghan 
doubt continue as computing power steadily increases. 

As with any numerical scheme however, SPH remains an approximation to real- 
ity, and as such reality checks are required in the form of standard tests. These act 
as calibration routines, to ensure the the results of any simulations are physically 
realistic, and can be relied upon. Many such tests exist, and there are far too many 
to do justice to here, but see for instance the astro code wiko, which has a number 
of cross-comparison tests with other codes, specifically aime d at d isc- like models. 



Pried (l2005i) the discussion 



As the code I use is a derivative of the one discussed in 
of numerical tests found here is particularly appropriate. A further suite of stan- 
dard tests including Sod shocks and Sed ov blasts among oth ers, used for both code 



verification and comparison, is given in ( iTasker et al.l . |2008| ). 



^http://"W"W"w-thcoric.physik.unizh.ch/astrosiin/codc/doku.php?id=homc:home 
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